
##clustercheckpint
#引用包
library(limma)
library(reshape2)
library(ggplot2)
library(ggpubr)

expFile="OV_TPM.txt"           #表达数据文件
ClusterFile="rs_tcga.txt"      #分型结果文件
geneFile="gene.txt"            #免疫检查点的基因列表文件
setwd("E:\\1.Ovarian_Treg\\1.picture\\10.checkpoint-diff")     #设置工作目录

#读取基因表达文件,并对数据进行处理
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

#读取基因列表文件，提取免疫检查点相关基因的表达量
gene=read.table(geneFile, header=F, sep="\t", check.names=F)
sameGene=intersect(row.names(data), as.vector(gene[,1]))
data=t(data[sameGene,])
data=log2(data+1)

#删除正常样品
group=sapply(strsplit(row.names(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[group==0,]
row.names(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", row.names(data))
data=avereps(data)

#合并数据
Cluster=read.table(ClusterFile, sep="\t", header=T, check.names=F, row.names=1)
sameSample=intersect(row.names(data), row.names(Cluster))
rt1=cbind(data[sameSample,], Cluster[sameSample,])
rt1$risk=ifelse(rt1$riskScore>median(rt1$riskScore),"high","low")
rt1=rt1[,c(sameGene, "risk")]

#提取显著差异的基因
sigGene=c()
for(i in colnames(rt1)[1:(ncol(rt1)-1)]){
  if(sd(rt1[,i])<0.001){next}
  if(length(levels(factor(rt1[,"risk"])))>2){
    test=kruskal.test(rt1[,i] ~ rt1[,"risk"])
  }else{
    test=wilcox.test(rt1[,i] ~ rt1[,"risk"])
  }
  pvalue=test$p.value
  if(pvalue<0.05){
    sigGene=c(sigGene, i)
  }
}
sigGene=c(sigGene, "risk")
rt1=rt1[,sigGene]

#把数据转换成ggplot2输入文件
rt1=melt(rt1,id.vars=c("risk"))
colnames(rt1)=c("risk","Gene","Expression")
#rt1=rt1[!rt1$Gene%in%c("BTNL2"),]
#绘制箱线图
bioCol=ggsci::pal_npg("nrc")(10)
#bioCol=bioCol[1:length(unique(rt1$risk))]
boxplot=ggboxplot(rt1, x="Gene", y="Expression", fill="risk",
                  xlab="",
                  ylab="Gene expression",
                  legend.title="risk",
                  width=0.8,
                  palette = c("#F26666","#212C5F"),outlier.shape = NA)+
  rotate_x_text(50)+
  stat_compare_means(aes(group=risk),
                     symnum.args=list(cutpoints=c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols=c("***", "**", "*", "")), 
                     label="p.signif",
                     label.y = 7.3)+
  ylim(0,7.5)+
  theme(panel.border = element_rect(fill=NA,
                                    color="black", 
                                    size=1, 
                                    linetype="solid"),
        axis.text.x = element_text(face = "italic"),
        axis.title = element_text(face = "italic"),
        legend.title = element_text(face ="italic" ),
        legend.text = element_text(face ="italic" ))#+coord_flip()

#输出图片
pdf(file="checkpoint.diff.pdf", width=12, height=4)
print(boxplot)
dev.off()

